%load_ext autoreload
%autoreload 2
%matplotlib inline
En este notebook se analiza en profundidad los resultados de las variables latentes obtenidas mediante el Autoencoder Convolucional. Para dar validez a estas variables se utlizarán como entrada a un algoritmo de clustering (KMeans), analizando cada grupo creado en base a su morfologÃa y su porcentaje de muestras de slides sanas y tumorales.
LibrerÃas
import yaml
import os
from tqdm import tqdm_notebook
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm
from matplotlib.patches import Patch
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from utils import plot_paired_imgs
from utils import read_images
from utils import imscatter
from model.cae import CAE
from keras.preprocessing.image import ImageDataGenerator
Configuración
with open('conf/user_conf.yaml', 'r') as f:
conf = yaml.load(f)
patches_path = os.path.join(conf['data_path'], 'slides', 'patches')
image_size = conf['wsi']['patch_size']
Se analizan conjuntamente los datos de train y test y que se obvervó que no habÃa ovefitting y que por tanto las variables generadas en ambos conjuntos serán similares.
train_df = pd.read_csv(os.path.join(conf['data_path'], 'train.csv'), sep='|')
test_df = pd.read_csv(os.path.join(conf['data_path'], 'test.csv'), sep='|')
data_df = pd.concat([train_df, test_df])
data_df.head(3)
Creación del lector de imágenes, como en este caso únicamente se va a a realizar la codificación el parámetro _classmode se pone a None.
imgen = ImageDataGenerator(rescale=1/255)
data_generator = imgen.flow_from_dataframe(data_df, directory=patches_path, x_col='filename',
target_size=(image_size, image_size), color_mode='rgb', class_mode=None, shuffle=False)
Carga CAE
Se carga uno de los modelos ya entrenados, en este caso se utilzará el modelo 1 que tenÃa una capa de 256 variables latentes.
model = 'model_1'
cae = CAE(path=os.path.join(conf['data_path'], 'models', model), load=True)
cae.latent_features
cae.filters
Predicción
Devuelve un array de dimensiones (72265, 256), es decir un vector de 256 elementos por cada una de las muestras.
data_enc = cae.encoder.predict_generator(data_generator,
steps=len(data_generator),
workers=8,
verbose=1)
data_enc.shape
256 variables es un número considerablemente alto como dimensión de entrada a un KMeans, este algoritmo al depender de distancias está muy afectado por lo que se conoce como la Maldición de la dimensión. Por este motivo se realiza primero una reducción a un espacio de 20 variables mediante PCA.
pca = PCA(n_components=30)
data_enc_pca = pca.fit_transform(data_enc)
Varianza explicada
Con las 30 primeras componentes de un PCA se explica el 60% de la varianza de los datos.
sum(pca.explained_variance_ratio_)
Aunque se puede obvervar que con las 3 primeras la varianza explicada ya es del 40%.
expl_variance = pd.DataFrame({'feature': range(1, pca.n_components + 1),
'explained_var': pca.explained_variance_ratio_,
'cumulative_explained_var': np.cumsum(pca.explained_variance_ratio_)})
expl_variance[['feature', 'explained_var', 'cumulative_explained_var']].head(10)
plt.figure(figsize=(12,5))
plt.bar(expl_variance['feature'] - 1, expl_variance['cumulative_explained_var'], align='center', alpha=0.7, color='navy')
plt.xlabel('Feature')
plt.xticks(range(pca.n_components), expl_variance['feature'])
plt.title('PCA Cummulative Explained Variance ')
plt.ylabel('Explained Variance (%)')
pass
A continuación se realiza un clustering con K-Means sobre las 30 componentes principales del PCA. Para ellos se calcula primero el número óptimo de clusters mediante lo que se conoce como el metodo del codo. Este método consiste en ejecutar el algoritmo sobre un conjunto limitado de Ks y calcular la incercia de cada uno.
La inercia es igual a la suma de las distancias cuadráticas de cada punto al centroide de su cluster, por tanto este valor siempre va a decrecer según aumenta K. Para elegir la K óptima se busca el punto en el que la cambia la pendiente de la curva Inercia/K, es decir, donde la disminución de la incercia empieza a ser menos pronunciada.
Selección de K (método del codo)
Ks = list(range(3,20))
inertias = []
for k in tqdm_notebook(Ks):
kmeans = KMeans(n_clusters=k)
kmeans.fit(data_enc_pca)
inertias.append(kmeans.inertia_)
Obvervando el gráfico no hay un codo muy claro, estarÃa entre 6 y 8. Puesto que en este caso una mayor separación puede ayudar a diferenciar mejor entre tipos de imágenes elige un K=8.
plt.figure(figsize=(10,6))
plt.plot(Ks, inertias, lw=2, marker='o')
plt.xlabel('K (number of cluters)')
plt.ylabel('Inertia')
plt.title('Elbow method for selecting optimal K')
Reentrenamiento con K óptima (K = 8)
K = 8
kmeans = KMeans(n_clusters=K)
data_cluster = kmeans.fit_predict(data_enc_pca)
data_df['cluster'] = data_cluster
Distribución de los clusters
El segmento 2 tienen un número muy bajo de muestras con respecto a los demás, esto quiere decir que tiene unas imágenes muy diferenciadas.
cluster_counts = data_df.groupby(['cluster']).size().reset_index().rename(columns={0: 'total'})
cluster_counts
Ejemplos de cada cluster
A continuación se dibujan unas pocas muestras de los patches en cada cluster, se aprecia que las imágenes se han agrupado por morfologÃa y color con bastante precisión. Destaca en especial el grupo 2, que contiene imágenes practicamente negras en su totalidad y hace pensar que son errores en la lectura.
Esta representación podrÃa ser etiquetada con ayuda de un médico especializado para poder hacer un análisis más exhaustivo.
for cluster in range(K):
print('\033[1m Cluster {}'.format(cluster))
sample_data = data_df[data_df['cluster'] == cluster].sample(20)
sample_imgs = read_images(sample_data['filename'], patches_path)
plot_sample_imgs(sample_imgs, n_rows=3, n_cols=6, size=2.5, color=True, shuffle=False)
plt.show()
Tipo de slide por cluster
Se analizan el ratio de tipos de slides en cada uno de los segmentos. Hay que recordar de nuevo que los patches etiquetados como Primary Tumor pueden ser patches sanos, pero no al revés.
cluster_types = data_df.groupby(['cluster', 'sample_type']).size().unstack().fillna(0)
total = cluster_types.sum(axis=1)
cluster_types['Primary Tumor'] = round(100 * cluster_types['Primary Tumor'] / total, 2)
cluster_types['Solid Tissue Normal'] = round(100 * cluster_types['Solid Tissue Normal'] / total, 2)
cluster_types
Estas distribuciones pueden orientar a identificar qué patches son realmente tumorales y cuales sanos, aquellos que pertenezcan a grupos con porcentajes altos de muestras sanas como los grupos 4 y 7 podrÃan ser idenficados como sanos. Los resultados aún asà no son concluyentes y requieren un análisis en mayor profundidad, por ejemplo se deberÃa analizar en cuántos clúster caen los patches de cada slide.
data_df.groupby(['cluster', 'sample_type']).size().unstack().plot(kind='bar', stacked=True, figsize=(12,7), cmap='coolwarm_r')
plt.title('Sample type distribution per cluster')
pass
Porcentaje medio de tipos de células por cluster
Un análisis similar al anterior se puede hacer calculando el porcentaje medio de los tipos de células de las slides de cada segmento. Se esperarÃa ver grupos con porcentajes muy diferenciados. Las conslusiones son parecidas a las anteriores, aunque los porcentajes de células no afectadas son mayores en grupos 4 y 7.
data_df.groupby(['cluster'])[['percent_normal_cells', 'percent_stromal_cells', 'percent_tumor_cells']].mean()\
.plot(kind='bar', stacked=True, figsize=(12,7), cmap='coolwarm')
plt.title('Average distribution of cell types in each cluster')
Representación de las muestras de cada cluster en 2D de PCA
Por último se representan la distribución de un conjunto de muestras en las dos primeras componentes principales del PCA y el clúster al que corresponden. Este dibujo ayuda a ver qué representa cada una de las componentes del PCA. Por ejemplo la primera componente, el eje X, parece identificar los colores oscuros y la segunda componente los espacios blancos; las imágenes con valores más altos en esta componente son opacas.
data_df['pca1'] = data_enc_pca[:,0]
data_df['pca2'] = data_enc_pca[:,1]
data_df['paths'] = data_df['filename'].map(lambda x: os.path.join(patches_path, x))
fig, ax = plt.subplots(figsize=(25,20))
cmap = cm.get_cmap('Set1')
n_samples = 30
legend_elements = []
for cluster in range(K):
if cluster in [2]:
continue
color = cmap(cluster / K)
samples = data_df[data_df['cluster'] == cluster].sample(n_samples)
imscatter(samples['pca1'], samples['pca2'], samples['paths'], zoom=0.4, ax=ax, color=color, lw=6)
legend_elements.append(Patch(facecolor=color, edgecolor=color, label='Cluster {}'.format(cluster)))
ax.legend(handles=legend_elements, loc='upper left', prop={'size': 15})
plt.axis('off')
plt.title('Representación de los clúster en 2D (PCA)', fontsize=20)
plt.show()